clear all

cd "$revdta"
set more off

set scheme s2mono
***************************Variable construction ***********************************************************************************************
capture program drop construct_vars_sp
program define construct_vars_sp
	xtset newid year
	sort newid
	merge m:1 newid using ".\intermediary\sp_temp_types"
	tab _merge
	drop if _merge!=3
	drop _merge
	sort newid year
	tab DS,gen(D)
	tab sp_DS, gen(spD)
	gen lag1D1=l1.D1
	gen lag1D2=l1.D2
	gen lag1D3=l1.D3
	gen lag1spD1=l1.spD1
	gen lag1spD2=l1.spD2

	/*pxy is the prob of moving from L=y to L=x */
	forvalues i=1(1)3{
	gen p`i'1= D1==1 if lag1D`i'==1 				
	gen p`i'2= D2==1 if lag1D`i'==1 
	gen p`i'3= D3==1 if lag1D`i'==1 
}
	forvalues i=1(1)2{
	gen sp_p`i'1=spD1==1 if married==1 & lag1spD`i'==1 
	gen sp_p`i'2=spD2==1 if married==1 & lag1spD`i'==1 
	}


end
******************************************************************************************************************************************



**************************Predict values for the <=62 part (to use in simulations)
foreach var in "" {
tokenize ${typename}_marginal 
*prodwage_marginalpost
foreach typetype in $type {
*type_prodpost
cd "$revdta"

use "Full_Panel_2b.dta", clear
keep if sample_health==1 & educ <= 2

bysort newid: egen firstage = min(age)
bysort newid: egen firsthealth = max(DS * (age==firstage))


local suffix

*keep if chd_maxhealth>=1 & chd_maxhealth<=3
*I lose most people before 1996 from child health, another 30% from restricting 
*to in-sample heads
*total: 2k households


*keep if marit==`relate'
*keep if married == 1 //I need the JOINT evolution of health...
gen disgroup_both =.
replace DS=DS+1
replace sp_DS = 1 if sp_DS ==2
replace sp_DS=sp_DS+1
replace sp_DS = . if married==0

replace disgroup_both=1 if DS==1 & sp_DS==1
replace disgroup_both=2 if DS==1 & sp_DS==2
replace disgroup_both=3 if DS==2 & sp_DS==1
replace disgroup_both=4 if DS==2 & sp_DS==2
replace disgroup_both=5 if DS==3 & sp_DS==1
replace disgroup_both=6 if DS==3 & sp_DS==2

replace DS=DS-1
replace sp_DS=sp_DS-1

construct_vars_sp
ren type_fix type_lp
ren `typetype' type_fix
levelsof type_fix, local(typevals)
***Probits by type and controlling for a full set of age dummies. Done for the baseline <=62 sample.
foreach j of local typevals{
	forvalues i=1(1)3{
	forvalues k=1(1)3{
	quietly gen  p`i'`k'_t`j'=.
	forvalues a = 23(1)62{
	quietly sum p`i'`k' if lag1D`i'==1 & type_fix==`j'&age==`a'
	quietly replace p`i'`k'_t`j'=r(mean) if lag1D`i'==1 & type_fix==`j'&age==`a'
					}
					}
					}
					
	forvalues i=1(1)2{
	forvalues k=1(1)2{
	quietly gen  sp_p`i'`k'_t`j'=.
	forvalues a = 23(1)62{
	quietly sum sp_p`i'`k' if lag1spD`i'==1 & type_fix==`j'&age==`a'
	quietly replace sp_p`i'`k'_t`j'=r(mean) if lag1spD`i'==1 & type_fix==`j'&age==`a'
					}
					}
					}

					
					}
	foreach x of varlist p??_t?{
	gen c_`x' =`x'
	}

	foreach x of varlist sp_p??_t?{
	gen c_`x' =`x'
	}
collapse (mean) p??_t? sp_p??_t? (count) c_p??_t? c_sp_p??_t?,by(age)
foreach x of varlist p??_t?{
replace `x' =. if c_`x'==0
}
foreach x of varlist sp_p??_t?{
replace `x' =. if c_`x'==0
}

drop c_p??_t? c_sp_p??_t?
keep if age>=23 & age<=62
sort age


foreach j of local typevals{
forvalues i=1(1)3{
egen rowtot`i'_t`j' = rowtotal(p`i'?_t`j'), missing
}
forvalues i=1(1)2{
egen sp_rowtot`i'_t`j' = rowtotal(sp_p`i'?_t`j'), missing
}
}

drop rowtot* sp_rowtot*

cd "$out"

save temp_2362,replace


**************************Predict values for the 62+ part (to use in simulations)

cd "$revdta"

use "Full_Panel_2b.dta", clear
keep if sample_select==1 & educ <= 2

bysort newid: egen firstage = min(age)
bysort newid: egen firsthealth = max(DS * (age==firstage))


local suffix

*I lose most people before 1996 from child health, another 30% from restricting 
*to in-sample heads
*total: 2k households

*keep if marit==`relate'
*keep if married == 1 //I need the JOINT evolution of health...
gen disgroup_both =.
replace DS=DS+1
replace sp_DS = 1 if sp_DS==2
replace sp_DS=. if married==0
replace sp_DS=sp_DS+1
replace disgroup_both=1 if DS==1 & sp_DS==1
replace disgroup_both=2 if DS==1 & sp_DS==2
replace disgroup_both=3 if DS==2 & sp_DS==1
replace disgroup_both=4 if DS==2 & sp_DS==2
replace disgroup_both=5 if DS==3 & sp_DS==1
replace disgroup_both=6 if DS==3 & sp_DS==2


construct_vars_sp
g agesel=age>=23 & age<=82
ren type_fix type_lp
ren `typetype' type_fix

preserve

gen part_full = lw!=. & hours >=1500 
gen part_part = lw !=. & hours < 1500 & hours >=500

gen sp_part_full = sp_lw !=. & sp_hours >=1500
gen sp_part_part = sp_lw != . & sp_hours < 1500 & sp_hours >=500

tab DS, gen(DSgroup)
*Marital-specific health risks
egen everwork = max(part_full), by(newid)
collapse (mean) p?? sp_p?? DSgroup? if insample_head==1 & type_fix==1,by(age married)
order age married p* sp_p*, alpha
sort age married
twoway connected p12 p13 age, by(married)
sum p11 if married==1
local mean1 = `r(mean)'
sum p11 if married==0
local mean0 = `r(mean)'
di "average diff: " `mean1'-`mean0'
restore

***Probits linear in age done for the 23-82 sample, assuming uniform growth by type after age 62
quietly{
	forvalues i=1(1)3{
capture probit p`i'1 age if lag1D`i'==1 & agesel==1
	if _rc == 0 predict q`i'1,pr
	if _rc != 0 {
	sum p`i'1 if lag1D`i'==1 & agesel==1 
	gen q`i'1 = r(mean) 
	}
capture probit p`i'2 age if lag1D`i'==1 & agesel==1
	if _rc == 0 predict q`i'2,pr
	if _rc != 0 {
	sum p`i'2 if lag1D`i'==1 & agesel==1 
	gen q`i'2 = r(mean) 
	}
capture probit p`i'3 age if lag1D`i'==1 & agesel==1
	if _rc == 0 predict q`i'3,pr
	if _rc != 0 {
	sum p`i'3 if lag1D`i'==1 & agesel==1 
	gen q`i'3 = r(mean) 
	}
}

forvalues i=1(1)2{
capture probit sp_p`i'1 age if lag1spD`i'==1 & agesel==1
	if _rc == 0 predict sp_q`i'1,pr
	if _rc != 0 {
	sum sp_p`i'1 if lag1spD`i'==1 & agesel==1 
	gen sp_q`i'1 = r(mean) 
	}
capture probit sp_p`i'2 age if lag1spD`i'==1 & agesel==1
	if _rc == 0 predict sp_q`i'2,pr
	if _rc != 0 {
	sum sp_p`i'2 if lag1spD`i'==1 & agesel==1 
	gen sp_q`i'2 = r(mean) 
	}
}
}					

collapse (mean) q?? sp_q?? ,by(age)
sort age

cd "$out"

merge age using temp_2362
drop _merge
keep if age>=23 & age<=82

****Now assumes transition probabilities for the 62+ grow using prob. linear in age (uniform growth by type) starting from flexible predictions by type at age 62
***Also deal with cases in which no flexible prediction at age 62 is available, so move back until flexible prediction is found
forvalues j=1(1)3	{
		forvalues k=1(1)3	{
			g g`j'`k'=q`j'`k'/q`j'`k'[_n-1]-1
							}
					}
					
	forvalues j=1(1)2	{
		forvalues k=1(1)2	{
			g sp_g`j'`k'=sp_q`j'`k'/sp_q`j'`k'[_n-1]-1
							}
					}


forvalues j=1(1)3	{
		forvalues k=1(1)3	{
			foreach t of local typevals{
				replace p`j'`k'_t`t'=p`j'`k'_t`t'[_n-1]*(1+g`j'`k') if age>62 | p`j'`k'_t`t'==.
								}
							}
					}

					
	forvalues j=1(1)2	{
		forvalues k=1(1)2	{
			foreach t of local typevals{
				replace sp_p`j'`k'_t`t'=sp_p`j'`k'_t`t'[_n-1]*(1+sp_g`j'`k') if age>62 | sp_p`j'`k'_t`t'==.
								}
							}
					}

***Finally, smooth (and interpolate to fill in empty age cells)
forvalues j=1(1)3	{
		forvalues k=1(1)3	{
			foreach t of local typevals{
				lowess p`j'`k'_t`t' age,gen(p`j'`k'_t`t'q) nog 
				ipolate p`j'`k'_t`t'q age,gen(p`j'`k'_t`t'p) e
				lab var p`j'`k'_t`t'p "Type `t'"
								}
							}
					}
					
forvalues j=1(1)2	{
		forvalues k=1(1)2	{
			foreach t of local typevals{
				lowess sp_p`j'`k'_t`t' age,gen(sp_p`j'`k'_t`t'q) nog 
				ipolate sp_p`j'`k'_t`t'q age,gen(sp_p`j'`k'_t`t'p) e
				lab var sp_p`j'`k'_t`t'p "Type `t'"
								}
							}
					}

*drop g* q* *q p11_t1-p33_t3
keep age p??_t?p p??_t? sp_p??_t?p sp_p??_t?
rename p??_t? rawp??_t?
rename sp_p??_t? rawsp_p??_t?
rename p*p p*
rename sp_p*p sp_p*
local reshapes
forvalues i=1(1)3{
forvalues j=1(1)3{
*collecting names for reshape
local reshapes `reshapes' p`i'`j'_t rawp`i'`j'_t
}
}

forvalues i=1(1)2{
forvalues j=1(1)2{
*collecting names for reshape
local reshapes `reshapes' sp_p`i'`j'_t rawsp_p`i'`j'_t
}
}

reshape long `reshapes', i(age) j(type)

forvalues i=1(1)3{
forvalues j=1(1)3{
*truncating negative interpolations at 0
replace p`i'`j'_t=0 if p`i'`j'_t < 0
replace rawp`i'`j'_t=0 if rawp`i'`j'_t < 0
}
}

forvalues i=1(1)2{
forvalues j=1(1)2{
*truncating negative interpolations at 0
replace sp_p`i'`j'_t=0 if sp_p`i'`j'_t < 0
replace rawsp_p`i'`j'_t=0 if rawsp_p`i'`j'_t < 0
}
}


rename p*_t p*
rename rawp*_t rawp*

rename sp_p*_t sp_p*
rename rawsp_p*_t rawsp_p*

reshape long p1 p2 p3 sp_p1 sp_p2 rawp1 rawp2 rawp3 rawsp_p1 rawsp_p2 , i(age type) j(to)
forvalues i=1(1)3{
bysort age type: egen p`i'tot = sum(p`i')
replace p`i'=p`i'/p`i'tot

bysort age type: egen rawp`i'tot = sum(rawp`i')
replace rawp`i'=rawp`i'/rawp`i'tot
}
drop p*tot rawp*tot

forvalues i=1(1)2{
bysort age type: egen sp_p`i'tot = sum(sp_p`i')
replace sp_p`i'=sp_p`i'/sp_p`i'tot

bysort age type: egen rawsp_p`i'tot = sum(rawsp_p`i')
replace rawsp_p`i'=rawsp_p`i'/rawsp_p`i'tot
}
drop sp_p*tot rawsp_p*tot


reshape long p rawp sp_p rawsp_p, i(age type to) j(from)

preserve
drop rawp rawsp_p
export delimited using "empiricalhealthprobs_`1'`suffix'", delimiter(",") replace
restore

gen bothhealth=from==1
gen spsick = to==2|to==4|to==6
gen headsick=to>=2
gen sevsick = to >= 3
collapse (sum) p rawp if bothhealth==1, by(age headsick type)

keep if age <= 62

twoway (line p age if headsick==1 & type == 1, lwidth(medthick)) (line p age if headsick==1 & type == 2, lwidth(medthick) lpattern(shortdash)) (line p age if headsick==1 & type == 3, lwidth(medthick) lpattern(longdash)), legend(rows(1) label(1 Low Type) label(2 Moderate Type) label(3 High Type) region(lstyle(none)) lwidth(none)) xtitle(Age, size(vlarge)) ytitle("Probability of Bad" "Health Shock for Ref. Person", size(vlarge)) ylabel(, labsize(vlarge))  xlabel(, labsize(vlarge))  plotregion(color(white)) graphregion(color(white)) 
/*
local combtypes
foreach t of local typevals {
local combtypes `combtypes' Type`t'
twoway (connected rawp age if spsick==0&headsick==0) (connected rawp age if spsick==0&headsick==1) (connected rawp age if spsick==1&headsick==0) (connected rawp age if spsick==1&headsick==1)if type==`t',  legend(label(1 Both Healthy) label(2 Head Sick) label(3 Spouse Sick) label(4 Both Sick)) xtitle(Age, size(vlarge)) ytitle(Probability, size(vlarge))  legend(off) yscale(range(0 1)) ylabel(, labsize(vlarge)) xlabel(, labsize(vlarge))  plotregion(color(white)) graphregion(color(white))
twoway (connected p age if spsick==0&headsick==0) (connected p age if spsick==0&headsick==1) (connected p age if spsick==1&headsick==0) (connected p age if spsick==1&headsick==1)if type==`t',  legend(label(1 Both Healthy) label(2 Head Sick) label(3 Spouse Sick) label(4 Both Sick)) xtitle(Age, size(vlarge)) ytitle(Probability, size(vlarge)) name(Type`t', replace) title(Type `t')  legend(off) yscale(range(0 1)) ylabel( , labsize(vlarge)) xlabel(, labsize(vlarge))   plotregion(color(white)) graphregion(color(white))
graph export "householdhealthshocks_type`t'_`1'.pdf", replace
}
grc1leg `combtypes', legendfrom(Type1) plotregion(color(white)) graphregion(color(white))
*/
graph export "householdhealthshocks_`1'`suffix'.pdf", replace

save "householdhealthshocks_`1'`suffix'",replace

macro shift
}
}
